Prictionless bead packs have macroscopic friction, but no dilatancy. 
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The statement of the title is shown by numerical simulation of homogeneously sheared assem- 
blies of frictionless, nearly rigid beads in the quasistatic limit. Results coincide for steady flows at 
constant shear rate 7 in the limit of small 7 and static approaches, in which packings are equili- 
brated under growing deviator stresses. The internal friction angle (f, equal to 5.76 ± 0.22 degrees 
in simple shear, is independent of average pressure P in the rigid limit and stems from the abil- 
ity of stable frictionless contact networks to form stress-induced anisotropic fabrics. No enduring 
strain localization is observed. Dissipation at the macroscopic level results from repeated network 
rearrangements, like the effective friction of a frictionless slider on a bumpy surface. Solid fraction 
$ remains equal to the random close packing value ~ 0.64 in slowly or statically sheared systems. 
Fluctuations of stresses and volume are observed to regress in the large system limit. Defining 
the inertial number as / = 'y y^m/JaP) , with m the grain mass and a its diameter, both internal 
friction coefficient n* — tSLiup and volume 1/$ increase as powers of I in the quasistatic limit of 
vanishing /, in which all mechanical properties are determined by contact network geometry. The 
microstructure of the sheared material is characterized with a suitable parametrization of the fabric 
tensor and measurements of coordination numbers. 

PACS numbers: 45.70.-n, 83.80.Hj, 81.40.Lm, 83.10.Rs 



I. INTRODUCTION 



Packings of particles appear in a variety of fields of 
condensed matter physics and material science, such as 
granular materials [11,11, Hkpowders [1| , or concentrated, 
non-colloidal suspensions a, @ . Such systems are macro- 
scopically disordered, and share many common features 
in their rheological behavior. One is a certain shear stress 
threshold, above which they roughly qualify as a fluid, 
and below which they might be regarded as solid. In 
assemblies of particles with purely repulsive force laws, 
interactions often do not introduce any stress scale, and 
the threshold only involves some ratio of stress compo- 
nents, whence a behavior often expressed as a friction 
law. Another basic property shared by many particulate 
systems is the existence of a specific value of the particle 
density, above which the material cannot flow. The vis- 
cosity of a dense suspension diverges as the solid fraction 
$ approaches some value $*, often regarded as identi- 
cal to the random close packing one, $rcp (^rcp — 0.64 
for identical spherical balls Q). Shearing and volumet- 
ric strains are coupled in granular media, which, once 
densely packed, cannot deform without expanding: this is 
the dilatancy property, first introduced by Reynolds 
Once the shear strain reaches a large enough value, gran- 
ular packs can continuously deform, like perfectly plastic 
materials, under constant stresses while keeping a con- 



stant solid fraction ^c'. this state of steady plastic flow 
does not depend on initial conditions and is known in 
soil mechanics as the critical state [lo| . Friction and di- 
latancy are coupled in granular materials by the stress- 
dilatancy relations, as proposed, e.g. by Rowe [Tol. [Tl|. 

It is tempting to try and identify simple, model systems 
apt to explore the microscopic origin of those broadly de- 
fined rheological features. To this end, discrete particle 
numerical simulation, for granular materials [H. 
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or suspensions |14| , has now become a widespread re- 
search tool. Thus friction laws in model granular as- 
semblies in steady shear flows, with some inertial ef- 
fects, were simulated [isl [T6j . and stress-dilatancy re- 
lations were tested [ITII . Many results were obtained 
on sphere packings [l8l. [l9|. which, long investigated in 
order to characterize their geometry are now 

studied with complete mechanical models. Thus it has 
been checked [2l|, [H, that the random close pack- 
ing state of monosized spheres is apparently uniquely de- 
fined if enduring agitation inducing traces of crystalline 
order is avoided in the assembling stage. The macro- 
scopic (or internal) friction coefficients, and their rela- 
tion to micromechanical parameters, including intcrgran- 
ular friction, have been evaluated from numerical simu- 
lations 

Despite recent advances, some open gaps and unan- 
swered questions can be pointed out in the literature. 
The accurate and detailed characterization of frictionless 
systems under isotropic loads [ll], in which static 
equilibrium states are studied, and few parameters are in- 
troduced, contrast with the more general investigations 
of the behavior of granular systems with intergranular 
friction [17, 18, 19], which are most often carried out by 
dynamical methods involving inertia effects, and involve 
quite a few additional parameters. In those latter stud- 
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ies, the limit of frictionless grains is not really treated 
with the desirable accuracy. Yet, frictionless packings, 
albeit reported to exhibit singular properties [Ij, HE, [la] , 
incorporate basic geometric effects that are common to 
suspensions and dry granular systems, even though they 
are supplemented by viscous flow effects in the first case, 
by intergranular friction and inertia in the second case. 

In order to clarify issues that have not been settled, the 
present paper is devoted to a numerical study of friction- 
less bead packings, subjected to homogeneous shear, and 
addresses the following questions. Can frictionless packs 
sustain shear stresses in static equilibrium states as well 
as in slow, steady flow, and do static and dynamic fric- 
tion coefficients coincide ? Do fluctuations on measured 
stresses or strain rates regress in the large system limit ? 
What can be said about characteristic densities $rcp ■, 
<i>c ? How do classical approaches of dilatancy [9l.l27j. 
and the way it couples to friction [13], apply in such a 
simple case ? 

The paper is organized in four main parts. Section |ll] 
describes the model material and the numerical simula- 
tion setup, specifying the boundary condition and initial 
states used in static and dynamic approaches. SectionlTTTl 
reports on the main results about the macroscopic behav- 
ior - macroscopic friction and dilatancy - and their de- 
pendence on the dimensionless control parameters identi- 
fied in Sectionini Section Hvl investigates the packing mi- 
crostructure and the force networks, in connection with 
macroscopic mechanical properties, with, in particular, 
a detailed characterization of anisotropy. Section IVlis a 
discussion. 



II. MODEL MATERIAL AND NUMERICAL 
EXPERIMENTS 

A. System and interactions 

We consider packings of equal-sized spherical beads of 
diameter a and mass to, enclosed in a cuboidal simulation 
cell. 

Beads interact in their contacts where only normal 
forces are transmitted, which are modeled as a sum 
of an elastic term and a viscous one, as in many numer- 
ical studies of granular systems (see e.g., Refs. [13, HI, 
m, H^). The elastic force is related to the normal 
deflection h of the contact by the Hertz law [30| . 

= f ^/^/^^/^ (1) 

where E is a. notation for E/{1 — v^), E is the Young 
modulus of the solid material the spherical grains are 
made of, and v its Poisson ratio. Eq. ^ attributes to 
contacts a variable spring constant Kn — dF^/Ah = 
E^h^'^/2. 

The viscous normal force opposes the relative normal 



velocity 5Vn = h of contacting beads, and is chosen as 

= C^2mKN5VN = CimEf'^iahf'HVN , (2) 

with a constant coefficient C. The same form of the vis- 
cous force was used in [HjIMI- Although ([2|) is devoid of a 
physical justification, some kind of dissipation is required 
(a granular material is not a conservative system), and 
consequently, the influence of C on the simulation results 
has to be carefully assessed. One attractive feature of the 
force law ^ and Q is the resulting velocity-independent 
coefficient of restitution e^v in binary collisions. Most 
simulations reported here were done with Q values such 
that Cat is close to zero. 

Particle rotations play no role and are ignored, as fric- 
tionless spherical objects behave like point particles in- 
teracting with central forces. 

The equations of motion for the particles, given by 
Newton's law, as in all molecular dynamics (MD) meth- 
ods, are to be numerically solved with standard time dis- 
cretization schemes fs^]. The time step used in the com- 
putations is a small fraction of the characteristic period 
of oscillations for the stiffest contact. 



B. Boundary conditions, stress and strain control 

We use different simulation procedures in which some 
strain, or strain rate, and stress components are exter- 
nally imposed to the system. 

In order to avoid wall effects and to determine eas- 
ily the intrinsic constitutive laws that apply in the large 
system limit, the simulation cell has periodic boundary 
conditions. The edges of the cell have lengths (iQ)i<Q<3 
along the three orthogonal axes of coordinates. Unlike 
the cell, the material may undergo some shear strain, 
imposed with the Lees-Edwards procedure [111 . Adding 
this possibility to the potential shrinking deformations 
along the three axes of coordinates, four independent 
strain components are considered in the different simu- 
lation steps and methods we are using in this work. The 
procedures defined below consist in choosing to fix some 
of them to zero or to a constant value while prescribing 
the values of stresses conjugate to the others. Table H] 
recapitulates those choices for the three different simula- 
tion procedures. 



1. Initial assembling process: procedure O 

In a preliminary step, the system is first prepared by 
isotropic compression of a loose "gas" of particles. The 
corresponding procedure, denoted as "O" (like "origin"), 
is the one applied in [23| to prepare isotropic packings. 
Global shear strain 7 is kept equal to zero, while the sys- 
tem shrinks along all three directions, until a mechanical 
equilibrium state is reached for which all three diagonal 
components Oaa of the Cauchy stress tensor (33l . |3J| are 
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equal to a set pressure value P. The system is deemed 
equilibrated when all forces compensate to zero, with a 
tolerance set to 10~'^ Pa^ on each particle, and each diag- 
onal stress component is equal to P with a relative error 
smaller than 10^^, while the kinetic energy per particle 
does not exceed lO^^Pa^. Those isotropic equilibrated 
configurations are the "random close packing states" , as 
studied in [Hill [H. 



Stress/strain control Procedure O Procedure D Procedure S 


CTll / Ll 


(711 =P 


constant Li 


CTll = P 


1^22 / L2 


0-22 = P 


C22 = P 


0-22 = P 


0"33 / L3 


0-33 = P 


constant Lz 


0-33 = P 


0"12 / 7 


7 = 


constant 7 


0\2 = T 



TABLE I: Choice of imposed stresses or strain rates in the 
three simulation procedures O, D, and S. 



2. Controlled shear rate: procedure D 

Initial configurations produced with method O may 
then be subjected to a simple shear deformation, in which 
a macroscopic motion along direction 1 is set up, with 
velocity gradient, on average, along direction 2 (by the 
Lees-Edwards procedure), while L3 and Li are fixed. L2 
is allowed to fluctuate in order to maintain (J22 equal to 
P on average (with very small fluctuations). The macro- 
scopic shear rate is denoted as 7. This defines procedure 
"D" (for dynamically sheared). It was implemented in 
a very similar way in [35]. One then records the time- 
averaged shear stress r = {au), as well as the sample vol- 
ume. It is important to note that Lees-Edwards bound- 
ary conditions arc fully compatible with either a linear 
velocity profile or very heterogeneous strain fields. With 
this procedure, shear strain 7 is set equal to the ratio of 
the offset along axis Xi of the neighbor copy of the sim- 
ulation cell in the X2 direction, to the length ^2- Conse- 
quently, due to fiuctuations in L2, the time derivative of 
7 is not strictly equal to 7 at all times. 



3. Static approach, controlled shear stress: procedure S 

In the limit of small 7, results of procedure D simu- 
lations should be comparable to static computations, in 
which the system equilibrates under an externally im- 
posed shear stress. To compare static and dynamic mea- 
surements (possible differences between "static" and "dy- 
namic" friction coefficients or threshold shear stresses 
in similar systems are mentioned in 16], and discussed 
in m, H^), we also implemented a completely stress- 
controlled, quasistatic procedure, denoted as "S" (for 
static). In procedure S, increasing values of shear stress 
T are stepwise applied, by increments 6t = 0.005 x P, to 
the initially isotropic configurations obtained with pro- 
cedure O, while the prescribed value of all three diagonal 
components aaa is the initial pressure P. 7, unlike in 
procedure D, is not constant. It satisfies a dynamical 
equation designed to impose a prescribed value t to (J12. 
For each value of r, one waits until a satisfactory equi- 
librium state is reached (with the same tolerance levels 
as in procedure O). The calculation is stopped when the 
packing does not equilibrate with the current value of r 
after 5 x 10^ MD time steps. The largest r value for 
which an equilibrium state was obtained is kept as an 
estimate of the shear stress threshold for onset of flow. 



I K C 

1 X IQ-'^ - 0.56 Ki = 3.9 X lO"*; k2 = 8.4 x 10^ 0.05 - 0.98 



TABLE II: Range of dimensionless parameters used in this 
study. 



C. Dimensional analysis, state parameters and 
geometric limit 

Assuming homogeneous steady states are observed in 
large enough samples, with shear rate 7 and normal stress 
P, then, by dimensional analysis [H, [l^, [2^ all dimen- 
sionless state variables, such as solid fraction 4> or average 
stress ratio {oyij <^ti) only depend on three dimensionless 
parameters. 

The first one, the inertial number, I — jy^^, charac- 
terizes the im por tance of inertial effects in dense granular 
flows lla . fig , [37^ and plays a central role for these sys- 
tems [3^, |3^- The quasistatic limit is the limit of / ^ 0, 
which we will systematically explore. 

The importance of contact deformation is character- 
ized by the second dimensionless parameter, a stiffness 
number which we deflne, like in [23], as k = (E/P)^^^. 
K is such that the typical contact deflection h under 
pressure P is proportional K^^a with a prefactor close 
to 1 [2^. In order to enable comparison of macro- 
scopic elastic properties with experimental results, we set 
E = 70 GPa and v = 0.3 (glass elastic constants). The 
pressure levels chosen, P = 10 kPa and P = 100 kPa, 
then respectively correspond to k = ki ~ 3.9 x 10^ and 
K = K2 — 8.4x10"^. These two values of k were reported 
to be large enough for the limit of rigid grains, i.e., of 
K —f +00, to be approached with good accuracy in the 
case of static packings ^2§\. 

Finally, the third dimensionless parameter is the level 
of viscous damping C, which appears in a viscous force 
and should not play a major role in the quasistatic limit. 

Table |TT] sums up the values of dimensionless control 
parameters used in the present numerical study. 

We should investigate the relations between global in- 
tensive variables, such as stresses, density, strain rate, in 
the limit of large samples, i.e., of ^ +00. It is ex- 
pected that for large enough samples the material state 
in shear flow will not depend on the specificities of bound- 
ary conditions, or on whether shear stress or strain rate is 
controlled. This requires the investigation of possible size 
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FIG. 1: (Color online) |cti2| (left axis, in black) and (T22 (right 
axis, in red) as functions of strain 7. Note that the left and 
right scales are different. Time series obtained with 7 = 3.2 x 
10"^ K = Ki, C = 0.98 and TV = 4000. 



effects and the study of the regression of fluctuations for 
global variables. Measured state variables should also be 
uniform in space - and thus one needs to check for possi- 
ble shear localization. If dimensionless variables such as 
stress ratios or density are well behaved in the triple limit 
of ^ CX3 (thermodynamic limit), / — > (quasistatic 
limit) and k, — * +00 (rigid limit), then the observed in- 
ner states and mechanical behavior of the packings only 
depend on their geometric properties - hence the name 
macroscopic geometric limit we adopted for such a situ- 
ation. One of the major goals of the present study is the 
investigation of material properties in this limit. 

Finally, as a practical application of the dimensional 
analysis of simulation parameters, let us note that the 
computational cost, expressed as a number of MD inte- 
gration steps needed to reach a given shear strain 7, is 
proportional to ^^/k/I. 



III. GLOBAL VARIABLES AND 
MACROSCOPIC BEHAVIOR 

Our global observations and measurements are re- 
ported in this section. Conditions for proper observations 
of the intrinsic behavior of the material subjected to pro- 
cedure D (shear-rate-controlled numerical experiments) 
are checked for in Section IIII A| in which various quali- 
tative aspects of the material state in shear flow are dis- 
cussed. Attention is then focused on macroscopic friction 
(Sec. IIIIBI) and dilatancy (Sec. IIII Cp properties of the 
material, which are more thoroughly and quantitatively 
investigated. Finally, the results obtained with proce- 
dure D at low inertial numbers are compared to those of 
the static approach, procedure S, in Section fill Dl Sec- 
tion |TTTE] discusses the essential results and their connec- 
tions with the literature on granular materials. 



A. Material state in slow shear flow: qualitative 
aspects 

With procedure D, we investigate steady states, and 
time series are collected for averaging. We are interested 
in intrinsic constitutive laws, as measured on averaging 
over the whole sample. It is therefore necessary to check 
for both invariance in time and homogeneity, in the sta- 
tistical sense. We should also assess the control of con- 
stant stress (722, and discuss the values of other stress 
components. 



1. Steady state flows and stress measurements 

Fig. [H displays the evolution of two components of the 
stress tensor, (T22 and CT12, with strain 7. It shows that 
(T22 is well controlled since it was requested to stay equal 
to E22 = 0.1 P in this mmrerical experiment. The evo- 
lution of stress 1712, from the initial, isotropically con- 
flned state, witnesses the existence of an initial transient, 
which has virtually ended at 7 = 0.1 in that case. The 
steady state part of the time series starts for values of 
strain 7 that depend on the inertial parameter, of order 
10^^ for the smallest / values, (~ 10~^), increasing typ- 
ically to about 0.5 for / = 10~^ and to several units for 
/ ~ 10~^. Unlike in dense systems with intergranular 
friction [i3i[iM[3j for which deviator stresses, starting 
from isotropically compressed initial states, go through a 
peak before approaching a plateau value at large strain, 
the shear stress in frictionless bead packs appears to grow 
monotonically, as a function of strain, toward its steady 
state value. Another notable feature of the shear stress 
as a function of time is the importance of fluctuations, 
which often exceed 30% of the mean value on the example 
of Fig. [2 in a sample of 4000 beads. A proper evaluation 
of average shear stresses thus requires careful statistical 
approaches and error estimates. 

As a practical criterion to detect the end of the ini- 
tial transient regime, we request that a small set of ba- 
sic measured quantities do not exhibit any visible trend. 
Specifically, shear stress ai2, volume fraction $ and coor- 
dination number z should all fluctuate about their mean 
value in a stationary manner, as well as the kinetic energy 
per particle, Scc, associated with velocity fluctuations. 
The latter is defined as 



SCc 



1 



N 



2N ^ 



(3) 



(5ec measures the instantaneous discrepancy between the 
actual flow generated by the Lees-Edwards boundary 
condition in the granular material and the affinc velocity 
fleld in a homogeneous continuum in shear flow. 

Unlike L2, lengths Li and L3 are constrained to remain 
constant in procedure D, so that an and txas may vary 
during the simulation. For / < 0.01, we observed that 
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time averages of an and (T33 differed from the initial hy- 
drostatic pressure P by less than 3%. This difference be- 
comes even smaller for smaller inertial numbers: for / = 
10-3, relative differences ((crii)/P) - 1 and ((0-33) /P) - 1 
respectively reduce to 1.0% and 2.2%. Those values de- 
crease down to 0.9% and 1.7% for / = lO^'', and to 
0.6% and 1.6% for / = 10"^. Although apparently not 
equal to zero, even in the quasistatic limit, those stress 
components are very small, and, consequently, will not 
be studied in the sequel. Sec. IIIIBi instead, focuses on 
accurate determinations of shear stress CT12. 

For a given number of particles, the relative fluctua- 
tions of the instantaneous value of (T12, $ and z (i.e., the 
ratio of their quadratic average to the mean value) seem 
to be independent of /. The average values of Sec, on 
the other hand, as compared to the kinetic energy of the 
macroscopic field, which is proportional to 7^, increases 
as / decreases. Fig. [2] is a plot of {6ec) /{ma'^j'^) ver- 
sus /, showing that this ratio approximately diverges as 
1// in the limit of / ^ 0. This agrees with measure- 
ments made in 2D simulations of shear flows: the same 
behavior is reported in Ref. [Tsj . and an interpretation 
was suggested, to which we shall return in Section IIII El 
These observations suggest that in the quasistatic limit 
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FIG. 2: Kinetic energy associated with velocity fluctuations, 
as defined in ([3]), normalized by ma^j'^, versus I, in simula- 
tions with 4000 beads, for k £ {k2, ki} and ( — 0.98. 

one has increasingly inhomogeneous instantaneous veloc- 
ity gradient fields, which wc now investigate. 



2. Instantaneous velocity profiles 

Instantaneous velocity profiles vi{x2) recorded at dif- 
ferent random times for different values of / are plotted 
in Fig.[3l Profiles Vi{x2) are obtained on averaging parti- 
cle velocities over slices cut alongside X2 in the simulation 
cell (particles overlapping slice boundaries contribute to 






FIG. 3: (Color online) Two velocity profiles at randomly cho- 
sen times, for 7 = 3.2 x 10"^ (top), 7 = 3.2 x 10"^ (middle), 
7 = 3.2 X lO"" (bottom). «: = ki, C = 0.98 and iV = 4000. 



several different averages). Inertial number / has an im- 
portant effect on the granular fiow. As shown in Fig. [3l 
instantaneous velocity profiles for / — 3.2 x 10"^ are lin- 
ear, whereas shear bands may appear for / — 3.2 x 10"^, 
as in the profile marked "L" (for "localized" ) on the bot- 
tom plot of Fig. [21 The transition between these two 
regimes seems to be gradual, with profiles in the middle 
part of Fig.m corresponding to / = 3.2 x lO^'', exhibiting 
somewhat intermediate features. 

Localization occurs here in the bulk of the material 
since the system is not enclosed between walls. Localiza- 
tion is thus an intrisic property of the studied material, 
which spontaneously appears for small values of / 40]. 

At first glance, it seems that the erratic behavior of 
the velocity profiles in the quasistatic limit may seriously 
jeopardize the interest of the results obtained on averag- 
ing over the whole simulation size and would demand 
specific analysis, distinguishing between material states 
within and outside shear bands. However, localization 
patterns are not persistent, and linear velocity profiles 
are recovered by time averaging, even in the J ^ limit, 
which means that on average, the flow is homogeneous. 
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FIG. 4: (Color online) Velocity profile after shear strain in- 
tervals 7 equal to 0.004 (red dotted curve), 0.02 (blue dashed 
curve) and 0.1 (black solid curve), following the instant corre- 
sponding to the localized profile marked L in Fig. [S] bottom 
graph. 

Figure m shows the gradual fading out of strain rate local- 
ization, after a strain interval of order 0.1. Shear bands 
thus randomly appear, move and disappear. Such a be- 
havior is witnessed by larger relative fluctuations of 5ec 
as / decreases. We did not carry out a detailed study of 
the lifetime and dynamics of nonpersistent shear bands, 
as the statistical homogeneity of the system in steady 
state shear justifies an analysis of global behavior based 
on averages over space and time. 



with the quadratic average of fluctuations of the observ- 
able quantity. In practice, due to intrinsic long-lasting 
correlations in the system, we observed that quite long 
runs were necessary. In some cases with / ~ 10~^, up to 
10^ simulation time steps (corresponding to a deforma- 
tion 7 ~ 4) were necessary for a correct evaluation of the 
uncertainty on /i*. 

In the present case, we could check that the time series 
of all observable quantities recorded in distinct samples 
differing only by their initial state were statistically iden- 
tical in steady state with high accuracy, as expected from 
critical state theory [lOl, 14^, S, Hi] . 

/X*, as estimated from time series in type D simula- 
tions, may depend on the three dimensionless numbers 
introduced in Sec. Ill CI as well as on the number of 
particles. This dependence is investigated in the follow- 
ing paragraphs. 

1. Effect of I 

Among the three dimensionless parameters governing 
the behavior of the system, the inertial number / has 
the strongest effect on ji* . Fig. [5] plots ^* as a function 
of /. It shows that /i* is an increasing function of /. 
This dependence of the macroscopic friction coefficient 
on the inertial number is similar to the ones reported in 
the literatiire. as obtained by both simulations and ex- 
periments [isl [l^ . IstI , although most published results 
pertain to granular systems with friction in the contacts. 
Here ji* approaches a finite nonzero value /ljq in the qua- 
sistatic limit of / ^ 0, despite the absence of friction 
at intergranular contacts. /Iq coincides with the internal 
friction coefficient of the material in its critical state. 

Note the accuracy of the displayed curve: statistical 
uncertainties measured with the blocking method are 
comprised between 10^'' and lO^'^ and are thus invisi- 
ble on the graph. 



B. Macroscopic friction coefficient 

For D simulations, the macroscopic friction coefficient, 
which we denote as is set equal to the time average 
- in the steady state - of the ratio of the shear stress 
ai2 to the normal stress (T22 (we use a convention where 
compressive stress components are positive) 



m2\ 
C22 



(4) 



The simulations produce raw data in the form of time 
series. The steady part of the time series is isolated as 
explained in Sec. IIII Al and /i* can then be easily com- 
puted. To estimate the statistical uncertainty on the 
measurement of averages over finite time series, we use 
the "blocking" (or "renormalization group" ) technique 
presented in [4l| . This yields error bars on measurements 
of averages in finite systems which should not be confused 



2. Effect of K 

Near the rigid limit, the macroscopic behavior should 
refiect the absence of stress scale in the contact law: 
stress ratios and derived quantities such as the macro- 
scopic friction coefficient are expected to be indepen- 
dent of the average stress. The friction coefficient hardly 
changes between the two values of k used in our simu- 
lations, indicating that the rigid limit k — s- cx) is accu- 
rately approached. Those simulations were carried out 
with C = 0.98 for 1.8 x 10""* < / < 5.6 x 10"\ and 
the relative variation on fi* is less than 2% throughout 
this range of inertial parameter on varying the stiffness 
parameter from k — ki to k — K2. Thus we decided to 
gather the values obtained for the macroscopic friction 
coefficient with k = K2 and k = ki in Fig. O because the 
uncertainty on the macroscopic, geometric limit of /z* to 
be estimated will eventually exceed this small difference. 
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FIG. 5: Macroscopic friction /i* vs. inertial number / for 
stiffness parameter k € {k2, ki}, damping parameter = 0.98 
and number of beads A'' = 4000. The solid line is Eq. Q with 
the parameters of Table HVl (no visible difference on using best 
parameters with or ^2). 



3. Effect of C 

The viscous damping term is indispensable in the 
model, as the only source of dissipation, but its mag- 
nitude should be irrelevant in the quasistatic limit. Con- 
sequently, the influence of C on our results had to be as- 
sessed and we performed simulations for different values 
of / with C, — 0.98 (this value corresponds to a restitu- 
tion coefficient cat = 3.3 x IQ-^), ( = 0.55 (e^ = 0.17), 
C = 0.25 [CN = 0.49) and ( = 0.05 (ejv = 0.87). Our 
results show that for / < 10""^, the maximal relative 
variation of Q on the macroscopic friction /x* is less than 
1%. Far from the quasistatic regime, the influence of Q is 
no more negligible: the relative variation of /i* is greater 
than 10% on changing ^ when / > 10""'^. 



500 0.1f69 0.3f00 0.2f88 0.6367 0.0022 0.6403 

3.2 X fO*^ 1372 O.ffOf 0.f907 0.f609 0.6380 O.OOfS 0.6408 

4000 0.f090 0.f245 0.f43f 0.6387 0.0008 0.6404 

500" 0.1432 1.263 0.8378 0.6738 0.0178 0.7148 

3.2 X 10"* 1372 0.1209 0.1689 0.1779 0.6365 0.0016 0.6390 

4000 0.1197 0.1002 0.1519 0.6368 0.0010 0.6388 

500 0.1473 0.2091 0.2275 0.6316 0.0027 0.6360 

3.2 X 10"^ 1372 0.1457 0.1293 0.1966 0.6322 0.0016 0.6346 

4000 0.1458 0.0764 0.1752 0.6323 0.0009 0.6338 

500 0.2112 0.2045 0.3317 0.6193 0.0025 0.6234 

3.2 X 10^2 1372 0.2123 0.1197 0.2802 0.6197 0.0015 0.6223 

4000 0.2125 0.0694 0.2517 0.6200 0.0009 0.6215 



TABLE III: Influence of sample size A'^ on macroscopic fric- 
tion /i* and volume fraction "!> for different values of inertial 
number I, with n = n\ and C, = 0.98. Superscripts ", 100" 
denote the average of the top percentile values in the steady 
state part of the time series. 

"This numerical experiment displays shear-induced ordering, a 
feature observed only for systems of N < 1000 beads (see Ap- 
pendix |X] for details). 



The effect of the sample size on the macroscopic fric- 
tion is unnoticeable for N = 1372 and N — 4000 since 
the difference between the friction coefficients pertaining 
to these two sizes is less than the statistical uncertainty 
marring the accuracy on /Lt*. However, the impact of N 
on fi* cannot be neglected in the quasistatic limit for 
a system of iV = 500 beads. These results show that 
some minimum number of beads, of order about 1000, is 
required to approach the thermodynamic limit with an 
acceptable accuracy. 

The data of Table lllll also witness the regression of 
fluctuations of stress ratio /i* in the steady state in the 
large system limit. The results are compatible with the 
classical form for the decrease of fluctuations of collective 
variables, viz. Afi/ii* oc iV-1/2. Specifically, for / = 
3.2 X 10~^, K = Ki and ( = 0.98, a fit of the data to the 
following form: 



Afi/fi* = (7.6±0.7)iV"i/2 



(5) 



4. Effect of N 

The influence of the sample size on the average of the 
apparent friction coefficient |(Ti2|/(T22, was investigated 
on comparing results for three different numbers of par- 
ticles: N = 500, N = 1372 and N = 4000. Results 
are listed in Table IIIII We also recorded the standard 
deviations, denoted as A/i, and the average of the top 
percentile of the instantaneous values, denoted as ^*'^''°. 
Let us recall that we are dealing here with the fluctua- 
tions of the time series, which differ from the statistical 
uncertainties on the average values. 



has good statistical admissibility. This result is shown in 
Fig. [S] in graphical form (two additional sizes N = 2048 
and N = 2916 were also simulated). Therefore, we expect 
the steady state stress-strain curves such as the ones of 
Fig. [1] however noisy for the sample sizes simulated, to 
become smooth in the large system limit. 

5. Approach to the macroscopic geometric limit 

According to the previous parametric study, the geo- 
metric limit can be confidently studied as the limit of 



8 




FIG. 6: A^/^* as a function of for J = 3.2 x 10"^ n = ki 
and C, — 0.98. The solid line equation is relation ((S)). 



K 


Mo 


a 


c 


Kl 


0.101 ±0.004 


0.38 ± 0.04 


0.40 ± 0.07 


K.2 


0.100 ±0.003 


0.39 ± 0.02 


0.42 ± 0.03 



TABLE IV: Best fit parameters for Eq. (|6]) and the data ob- 
tained with — 4000, = 0.98 for k = ki and k = K2. 



/ ^ on samples of 4000 beads with k > K2 and 
C = 0.98. 

/i* should be a function of the sole inertial number in 
very good approximation for sufficiently small values of 
/. In the absence of any scale, we tried to fit the data by 
a power law function (see Fig. [5]) of the form 

fi* = f,* + cr (6) 

As stated above, this fit is not expected to be accurate 
for large / values and wc therefore restricted ourselves to 
fit the data points with / < 0.01. Parameters /ig (the 
geometric macroscopic friction coefficient), a and c were 
separately estimated for n = ki and k = K2 (keeping C = 
0.98 and N = 4000) and results are shown in Table HVl 

The value of the geometric macroscopic friction angle 
ifiQ corresponding to is (for k = ki) 

= 5.76° ± 0.22° (7) 

Quite similar values are also reported with two- 
dimensional packings of polydisperse disks by Taboada et 
al. [17], whose estimate of the macroscopic friction angle 
lies between 4° and 7° for frictionless grains, and by Da 
Cruz et al. p^ . who obtained /i^ ~ 0.1 in shear fiow 
simulations for small / parameters. Hatano [IQ] recently 
performed 3D numerical simulations on polydisperse as- 
semblies of about 10000 spherical beads, for different in- 
tergranular friction coefficients fj,. The reported value 
of the macroscopic friction coefficient in the quasistatic 
limit is 0.06 for /i = 0, apparently lower than our re- 
sult. It should be recalled however that Hatano's work 



was motivated by applications to granular materials un- 
der high confining stresses within geological fault zones, 
and that consequently simulations were carried out with 
lower stiffness levels (k = 1840, 136, 84 and 42) than in 
the present study. Moreover, the lower range of / pa- 
rameters, below 3 X lO"**, was only investigated with the 
lower K values. Hatano used the same form as Eq. © to 
fit his data, but his estimate a = 0.28 ± 0.05 differs from 
ours (see Table HV)). Although some effect of the poly- 
dispersity is possible, we also attribute this discrepancy 
to some non-negligible influence of k in Hatano's simu- 
lations [3]. Only the simulations with k = 1840 in [T^ 
can be expected to approach the rigid limit accurately. 
For this stiffness level, Hatano's data points are available 
for / > 3 X 10~^ and are in very good agreement with 
ours (e.g., n* ~ 0.17 for / = 0.01). 



C. Dilatancy and steady-state density 

Dilatancy under shear is a basic property of granu- 
lar materials in quasistatic deformation [l^, [HI, [13] j 
when dense samples are subjected to increasing deviator 
stresses. The steady-state density is mainly sensitive to 
/ if K is large enough [l^. The small / behavior of fric- 
tionless bead assemblies is investigated here with greater 
accuracy than in previous studies. 

We could check that, just like the macroscopic friction 
coefficient, the steady state time average of the volume 
fraction, $ = ($(t))t, is an intrinsic property of the ma- 
terial, independent of its initial preparation. Next, we 
investigate its dependence on the three dimensionless pa- 
rameters /, K and ^ and on the number of particles N. 



1. Effect of I , K and 

Once again, numerical experiments demonstrate that 
among the three dimensionless numbers governing the 
behavior of the system, the inertial number / has the 
most important effect on $ = {^{t))t. Fig. [7] shows the 
influence of / on $. We observe that $ decreases for 
increasing /, as previously reported [H, [l3|. It starts 
from a value <I>o — 0.64 in the quasistatic limit and the 
system expands as / increases. Statistical uncertainties 
on <i> measured thanks to the blocking method are com- 
prised between and lO^** and are thus invisible on 
the figure. $0 is very close to $rcp [HI HI], which co- 
incides (up to small corrections due to the finite con- 
tact stiffness) with the initial volume fraction cE>'*'°, right 
after the samples are assembled with procedure O. For 
K = K2 and iV = 4000 we have = 0.6382 ± 0.0011, 
and = 0.6369 ± 0.0009 for k = ki (averages and 
standard deviations are evaluated on five samples). The 
system studied thus appears to be devoid of dilatancy 
in the quasistatic limit. Whether $0 should be regarded 
as equal to ~ $rcp at the macroscopic level will be 
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0.64 



0.62 



0.6 



0.58 




Q I I I I I I I 

10"^ 10"" 10"^ 10"^ 10"' 10° 

I 

FIG. 7: (Color online) Volume fraction $ as a function of 
inertial number I (for ( = 0.98, N = 4000), for both stiff- 
ness levels K = Ki (blue squares) and k = K2 (red triangles) 
The solid lines are given by Eq. with the parameters of 
Table IVl 



uL 



uL 



discussed after the possible influence of N on the average 
densities is investigated. 

Stiffness parameter k typically induces a relative in- 
crease of the volume fraction of roughly 0.1% when it 
changes from k = K2 to k = ki, whatever the value of / 
- a small effect, yet distinguishable from statistical un- 
certainties. Such a density increase is of course expected, 
as larger contact deflections due to larger stresses or a 
softer material decrease the sample volume. Simulations 
with C = 0.98 (e^v = 3.3 x IQ-^), C = 0.55 {cn = 0.17), 
C = 0.25 (ejv = 0.49) and C = 0.05 (ejv = 0.87) for a wide 
range of inertial numbers have also been run. The influ- 
ence of C on $ is important for large /: for / > 0.1, the 
relative variation of with C can reach 30%. However, 
this effect, as expected, gradually vanishes as the qua- 
sistatic limit is approached, and for / < 0.01 the relative 
variation of $ with C is less than 0.1%. 



2. Effect of N 



According to Table IIIIl $ very slightly varies with 
the number N of particles, like in static, isotropic sys- 
tems [21I, [2^. The following fit, based on the measure- 
ments for the smallest available value of /, i.e., I = Ii = 
3.2 X 10^^ for K = Ki, may be used: 



$(k ^ Ki,I ^ h,N) = $1 - kiN-^/'^ 




FIG. 8: A$/$ as a function of A'^ for the same time series as 
in Fig. H fitted with Eq. ^ (solid line). 



with the following parameters: 



$1 = 0.6398 ±0.0002 
ki = 0.070 ±0.008 



(9) 



As with the macroscopic friction coefficient jj,* , we could 
check for the regression of fluctuations of the volume frac- 
tion for increasing N. For the same set simulations with 
/ = 3.2 X 10"^ K = Ki and C = 0.98 as in Sec. IIIIB4| 
we observe that the decrease of density fluctuations with 
increasing N can be fitted by the following relation: 



A 

N' 



with A = 0.051 ±0.011 



(10) 



as shown graphically in Fig. [8l whence a well-defined 
large system limit for $. 



3. Approach to the macroscopic geometric limit 

Volume fraction $ can therefore be modeled as a func- 
tion of / near the quasistatic limit, say for / < 0.01, with 
small corrections to account for the infiuence of N and 
K. It can be regarded as independent of C (at least for 
/ < 0.01). The fit form used is 



er 



(11) 



For N = 4000 and ( = 0.98, the best fit values of the 
parameters of Eq. (jlip are given in Table |Vl 

To evaluate the macroscopic value $0 ™ double 
limit of / ^ and A'' ±00, it is reasonable to assume 
that the small corrections to $g that result from the fi- 
nite value of N and from the nonvanishing value of / are 
additive. The use of Eqs. (0)-® to evaluate the finite A^ 
correction to the value $0 of the quasistatic density, as 
obtained on fitting Eq. to the results with A^ = 4000, 
yields, for k = ki: 



(8) 



0.6410 ±0.0005 



(12) 
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K <l>o V e 





0.6398 ± 2 


10-* 


0.39 ±0.01 


0.1786 ±8. 10"* 


K2 


0.6405 ± 2 


10-* 


0.42 ±0.02 


0.2038 ±3 10"* 


TABLE V: Best fit paramet 


ers for Eq. (Illf) and tiie data ob- 


tained with — 


4000, C = 


= 0.98 for K = 


/^i cLnd hi = hv2 ■ 


N 




Sn 


(At"'"*) 




256 




4 


0.246 


0.022 


500 




4 


0.210 


0.007 


1372 




6 


0.169 


0.004 


2048 




6 


0.154 


0.004 


2916 




6 


0.145 


0.007 


4000 




10 


0.136 


0.007 


8788 




6 


0.122 


0.005 



TABLE VI: Average (^''*'"> and standard deviation A/i='"' of 
tlie static friction coefficient obtained in S-type simulations, 
over 5'jv samples of N grains for different A'^. Data corre- 
sponding to both values of k (with Sjv /2 samples each) are 
aggregated. 
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0.01 0.02 0.03 

FIG. 9: (Color online) Size dependence of {^j,"*''^*). N denotes 
the number of particles in the system. The solid line is the 
fit of Eqs. p4|l - (|15p . Crosses are the top percentile values 
extracted from the time series of |(ti2|/(722 obtained in proce- 
dure D, as listed in Table Hill The hashed region represents 
the estimate, from D simulations, of fig with its error bar (Ta- 
ble |lVll. The (blue) dot with an error bar on the left axis is 
the static estimate, /i^?*. 



The increases of $, from its value in the rigid limit, due 
to the finite stiffness is of order (see [23l Eq. 31]) 
and is smaller than the statistical uncertainty in (|12p . 
The value of $q given in is thus our best estimate, 
from D-type simulations, of the solid fraction of sheared 
sphere packings in the macroscopic geometric limit. 



D. Static behavior 



We now compare the results of Sections HITB] and HITC] 

for steady shear-rate-controlled simulations (procedure 
D) with those obtained through static shear numerical 
experiments (procedure S). 



of frictionless disks: in the limit of vanishing shear rates, 
the shear stress reaches its large system limit with only 
several hundreds of beads, whereas the minimum shear 
stress required to maintain a long lasting steady shear 
flow exceeds the previous one and is more sensitive to N. 

Fig. [9] shows the influence of system size N on fj,^*^^^ 
(discarding the smallest N values). The data are cor- 
rectly fitted by the following relation 



with 



/i^^' = 0.091 ±0.009 



2.87 ±0.32 



(14) 



(15) 



1. Friction coefficient 

The static macroscopic friction coefficient is defined in 
procedure S as 



(13) 



P 



system has been able to sustain in mechanical equilib- 
rium, and P the confining pressure. Static microscopic 
friction coefficients for different sample sizes are displayed 
in Table IVII Values of ff^^^ are larger than the dynami- 
cal value fiQ — 0.100 ±0.004 obtained in D simulations in 



the quasistatic limit. As shown by Tab. IVH ii is size 
dependent, unlike /ip (for N > 1000). Analogous obser- 
vations were reported in [29j for two-dimensional systems 



The related angle of friction is (p^^* = 5.20° ±0.52°. This 
is consistent with the equality, in the thermodynamic 
limit {N oo), of the dynamical and static macroscopic 
friction coefhcients (see Eq. [7]). The influence of k is 
very small and negligible in comparison to the effect of 
the system size, and we therefore averaged over systems 
with both stiffness levels ki and K2- 

With the smallest system size simulated, N — 256, 
we observed that some of the samples, once submitted 
to shear stresses, acquired a strongly ordered, crystalline 
structure, to be discussed in Appendix \X[ 



2. Density 

Static shear simulations with procedure S support the 
observation made in Sec. IIIICl that the frictionless model 
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0.641 




|t|/P 

FIG. 10: Variation of the volume fraction $ with the static 
shear \t\/P imposed to five different samples of 4000 beads 
with K = K2. Each curve stops at a given value of \r\/P — 
|Tmax|/-P: this is the greatest value for which the packing has 
managed to reach mechanical equilibrium. 



fractions in quasistatic shear are identical, within statisti- 
cal uncertainties. Disregarding the very small correction 
due to the finite value of ki (equal to about 1.1 x 10^'' on 
applying the formula given in f23', Eq. 31]), this means 
that, just like for fi* , the values of solid fraction $ in 
the macroscopic, geometric limit coincide in strain rate 
controlled and in shear stress controlled approaches. 

As to the value of the solid fraction in the initial 
isotropic state, a similar evaluation of size effects yields 
(using the samples of Table IVII with k = ki and N > 
500): 

with 

= 0.6397 ±0.0008 , . 

ko = 0.15 ±0.03. 

As announced, this is the random close packing value [2l|, 
[23l |. Results ([19]), ^7}i and ((T2]) are compatible, and we 
thus conclude that the system is devoid of dilatancy un- 
der shear in the macroscopic geometric limit. 

E. Discussion 



material studied is devoid of dilatancy in the quasistatic 
limit. As shown in Fig. [101 which represents $ as a func- 
tion of the macroscopic stress ratio t/P imposed to the 
material in different samples with N = 4000, the volume 
fraction hardly evolves with the stress deviator as it is 
increased towards its maximum value. However, what- 
ever the initial state of the system, it experiences a slight 
compaction at the beginning of the shear and a small 
decompaction near the failure limit, but we have no con- 
vincing explanation for this phenomena. The evolution 
of $ is somewhat erratic (as in previous studies on 2D 
rigid, frictionless disk assemblies \24 . i25j ) and the den- 
sity change between the isotropic initial state and the one 
supporting the maximum shear stress is equal to zero, 
within statistical uncertainties. Similarly to the values 
of <1> measured in D simulations, solid fraction in 
static packings under maximum shear stress is slightly 
dependent on sample size, with a negative finite-size cor- 
rection to the macroscopic value. On fitting a variation 
proportional to N~-^/^ one gets, for k — ki, 

with 

= 0.6403 ±0.0004 
°° 17 

k = 0.125 ± 0.026. ^ ^ 

$ values for k — K2 are slightly larger, by about 10^''. 

Comparing this estimate of with the result for <i>g 
given in (|12p. we conclude that static and dynamic solid 



We briefly review and comment here the essential re- 
sults on the macroscopic behavior of the material under 
simple shear, and compare them to other available results 
on similar systems. 

1. Internal friction and the macroscopic geometric limit 

Whether assemblies of frictionless grains have a well- 
defined, finite internal friction coefficient has sometimes 
appeared as a debatable issue, although some previously 
cited works [isl . [l6l . [45| relying on numerical simulations 
of slow shear flows in steady state agree with our positive 
conclusion. A proper evaluation of /Iq in the macroscopic 
geometric limit requires more care than corresponding 
measurements in granular assemblies with friction. This 
is due to the importance of fluctuations, as apparent on 
Fig. [TJ In D-type simulations, it is also necessary to ex- 
plore a range of very small inertial numbers to accurately 
evaluate the quasistatic friction coefficient, as apparent 
in Fig. [51 As an example, for / = 5.6 x 10"^, quite a 
small value, the macroscopic friction coefficient already 
exceeds its quasistatic limit by 25%. 

Our estimate for /Iq is confirmed by static simulations, 
once the results are suitably extrapolated to the limit of 
large systems. One may understand this size effect on 
S results as follows. The friction coefficient evaluated in 
D simulations is an average over time series with large 
fluctuations. However, the system remains close to me- 
chanical equilibrium at any time. Assuming it is possi- 
ble to find an equilibrated configuration very close to all 
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dynamically explored states, however large the instanta- 
neous value of the shear stress, the S procedure would 
be able to find statically supported shear stress values as 
large as the maximum of ai2 in D time scries. Although 
clearly oversimplifying the evolution of the system in con- 
figuration space, this explanation appears to be correct at 
least on correlating A'^-dependent maximum static shear 
stress levels to fluctuations in slow shear flows: the N 
dependence of /z'^*'^*, as plotted in Fig.[9l is paralleled by 
that of the typical largest values of ji* (top percentile) in 
D simulations. 

Static and dynamic values of shear stress thresholds 
for flow are also observed to coincide in the flxed density 
simulation results of Xu and O'Hern [2^, obtained on 2D 
packings of frictionless disks, with a similar excess of the 
static value that vanishes as TV ^ cxd. 

When non negligible inertial effects are present, we ob- 
serve that the increase of internal friction with inertial 
number / is the dominant feature of the material behav- 
ior (the effect of stiffness level k is smaller by orders of 
magnitude), in qualitative agreement with many other 
results on frictional and frictionless grains ^Ji, tl6|| . 



2. Absence of dilatancy 

Our results also agree in static shear and in steady 
state constant shear rate flows for the average volume 
fraction which stays equal to its value in the initial, 
isotropically confined configuration in the macroscopic 
geometric limit. Our data show that, within statistical 
uncertainties (i.e., about 5 x 10"^) the critical value of 
$ is equal to <1>rcp in packings of frictionless spherical 
beads. 

The material studied is thus devoid of dilatancy. Inter- 
estingly, this contradicts the simple pictures of the origins 
of dilatancy which have been proposed since the intro- 
duction of this property by Reynolds [9*1, based on the 
distortion of simple assemblies of a small number of con- 
tacting spheres (like, e.g., a regular tetrahedron) |27]. 

The absence of dilatancy in the quasistatic limit is also 
at odds with the classical ideas on the relation between 
dilatancy and internal friction, according to which macro- 
scopic friction stems from two microscopic origins, inter- 
granular friction and dilatancy, with an additive combi- 
nation of relevant angles P,[23. Ref. [i3 adds another 
component ipo to macroscopic friction, due to intergranu- 
lar collisions as a source of dissipation, and therefore ac- 
counts for the internal friction of frictionless grains. Thus 
(/?o is the internal friction angle that we measure in the ge- 
ometric limit. Ref [l3|, although only incidentally dealing 
with frictionless materials, nevertheless appears to pre- 
dict a positive dilatancy in that case, which our results do 
not confirm. Similarly, a recent study published by Kruyt 
and Rothenburg [ij, which also deals with 2D disk as- 
semblies, predicts a non-vanishing dilatancy when inter- 
granular friction coefficient approaches zero. Ref. [i^ . 
similarly to Ref. [l3|, discusses stress-dilatancy relations. 



and finds a linear variation of the dilatancy ratio with the 
difference between peak and steady-state macroscopic 
friction. In contradiction with our data, it attributes 
a positive value to both quantities as the friction coef- 
ficient approaches 0, while its estimate for /ip is signif- 
icantly larg er than our (3D) one, or than the (2D) one 
of Ref. (171. (Note that the maximum deviator to mean 
stress ratio, as defined in J^, 44], is smtf ~ tamf). The 
frictionless case was not directly simulated in this work. 
Some rapid variations of macroscopic friction and dila- 
tancy angles near the singular limit of /xq ^ might be 
overlooked. 

In our simulations, instantaneous fluctuating shear 
stress and volume fraction, however, appear to be cor- 
related, suggesting some stress-dilatancy coupling at the 
level of short-lived, transient and rearranging structures, 
which disappears on taking time averages. 



3. A toy model 

Since some of our results on the macroscopic behavior 
of frictionless bead packs might seem counter-intuitive, 
we designed a simple model in which similar basic in- 
gredients (geometric constraints defining isolated equi- 
librium positions, inertia, viscous dissipation) produce 
an analogous behavior in a suitably defined "macroscopic 
geometric limit" . In both cases, the microscopic motion 
is a succession of arrested dynamical phases, alternating 
with approaches to transient equilibria. The toy model 
simply provides suggestive analogies, it should not be re- 
garded as a real physical explanation for the macroscopic 
behavior of the granular system. 

We consider a single object of mass Af , subject to its 
weight W ^ pushed along a rough horizontal surface. For 
simplicity the model is two-dimensional, with only one 
horizontal coordinate, xi, and the surface profile h{xi), 
along vertical coordinate X2 , is periodic with wavelength 
A, as depicted in Fig. [TT] The mobile object is driven 
either by a constant horizontal force F, or by a piston 
with constant horizontal velocity V. Both contacts are 
rigid and unilateral, so that the mobile object might move 
faster than the piston if accelerated downhill by gravity, 
or occasionally take off from the surface. Force F is the 
analog of shear stress ai2 in the granular material, and W 
that of (722, while horizontal and vertical displacements 
respectively correspond to shear strain and volume in- 
crease. A viscous force opposes the tangential motion 
along the surface, so that for F = the slider stabilizes 
at some local minimum of profile h{xi), like point O on 
the figure. Such minima are analogous to equilibrium 
states under isotropic pressure. 

Let us first discuss the static experiment, in which the 
mobile object, initially in equilibrium in O under F — 0, 
is subjected to a growing horizontal force F. It equi- 
librates where the tangent direction to the substrate is 
orthogonal to the applied force, ^ = F/W. It has first 
to move upwards, hence some dilatancy. The maximum 



13 




FIG. 11: The model of the shder on a rough surface, (a) Case 
of a sinusoidal profile. Forces at point S are drawn as vectors, 
(b) Profile for which fio ~ 

value of F/W is the static effective friction coefficient 
fis — tan (p of the object on the surface, equal to the 
maximum slope of profile h(xi). It is reached at point S 
on the figure. Effective static friction angle ip is the max- 
imum angle between the reaction of the substrate, force 
R on Fig. [TI] and the vertical direction. As the qua- 
sistatic motion from O to S follows the surface, dilatancy 
tan ■)/), defined as the ratio of vertical to horizontal coordi- 
nates of the velocity (corresponding to ratio ^22/7 in the 
sheared granular material), is also identical to the maxi- 
mum profile slope. Dilatancy and friction angles coincide: 
tjj = (p. If a nonzero friction coefhcient /io = tan ipQ is in- 
troduced in the contact between the mobile object and 
the substrate, then reaction R at point S (see Fig. [TT|) 
may form an angle ipo with normal direction (Sn), so 
that the effective static friction angle is f = ^0 +_0 - a 
classical form of the stress- dilatancy relation [iCt Il7|. 

In the velocity-controlled case, the dynamic friction co- 
efficient is conveniently evaluated from the dissipation of 
energy. In the limit of small velocity, the mobile object 
pulls ahead of the velocity-controlled driving piston at 
each maximum of h{xi). Its subsequent downhill sliding 
is accelerated by gravity, but it is prevented by viscous 
dissipation to pass the next maximum, and ends up at 
the bottom of the valley, where it is later picked up by the 
slow piston, to be pushed up the next ascending slope. 
In this scenario the dissipated energy per wavelength A 
is the potential energy loss HW in a fall over height H . 
Hence an effective friction coefficient hd = H/X. This 
result is, remarkably, independent of the viscous damp- 
ing coefficient, just like the macroscopic friction of the 
granular material in slow shear flow. As the properties 
of the system only depend then on substrate geometry. 



the limit of slow imposed velocity is the geometric limit. 

The macroscopic limit can be defined as X/L — > 0, 
where L is the length scale on which the effective prop- 
erties of the slider are studied. Consequently, its vertical 
motion, on the scale ~ A of microscopic asperities of 
the surface, becomes irrelevant, and one observes effec- 
tive macroscopic friction without dilatancy. Models for 
dilatancy [13] apparently focus on microscopic phases of 
the motion in which the slider rises up the slope, but 
ignore the equally important ones in which it falls down. 

Hd is the average slope of the ascending part of the pro- 
file, multiplied by the fraction of length for which h{xi) 
is increasing. It is in general smaller than /ig, which 
is the maximum slope. Thus, for a sinusoidal profile, 
h{xi) — H /2sin{2iTXi/ X), as represented on Fig. [TlTa). 
one has = t^H/X, while /id = H/X. In order for both 
friction coefhcients to coincide, function h{xi) should be 
as shown on Fig. Illf b). with ascending parts of constant 
slope, followed by vertical drops. 

The particular profile shape of Fig. [TlTb) can be ar- 
gued to be appropriate for the analogy with the bulk 
material. As long as the contact with the substrate is 
maintained, the configuration might be an equilibrium 
position for some (possibly negative) value of F . In the 
analogy with the granular material, the contact network 
might balance the external load for some value of the 
applied stress components. The free fall, on the other 
hand, is the analog of a network rearrangement, during 
which applied loads cannot be supported because of the 
missing contacts. In the granular material (as explicitly 
shown in [131) intervals of stress components for which 
a given contact network is stable shrink to zero in the 
large system limit. This corresponds for the toy model 
to a constant slope of the rising parts of profile h{xi) 
(defining a unique possible value of F in equilibrium) . 

Finally, the velocity-controlled sliding of the object 
on the profile of Fig. Illf b) also provides an interpre- 
tation of inertial number I (37| . and of the behavior 
of the kinetic energy [l5| . The motion involves two 
characteristic times: the duration of the rising phase, 
in which the object is in contact with the piston and 
moves with horizontal velocity V, ti = X/V; and the 
duration of the free fall, T2 cx y/ (MH/W). Their ra- 
tio, T2/T1 oc (V/X)^ (MH/W), is the analog of number 
/, as readily checked on replacing distance H by some 
length of order a (a typical interstice between neighbor- 
ing grains to be closed for a new contact to appear), 
V/X by 7, and force W by Pa^, which is the order of 
magnitude of unbalanced forces on the grains in the dy- 
namical phases of motion. The free fall phases of the mo- 
tion explain why the kinetic energy is, on average, much 
larger than the scale MV'^ associated with the macro- 
scopic motion. More precisely, the time average Scc of 
the kinetic energy associated with velocity fluctuations 
is of order HW{t2/ti) (for T2 ^ ti), whence (discarding 
constant factor {H/XY) the behavior shown in Fig. [21 
See cx {MV^)/I. 
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IV. MICROSTRUCTURE AND FORCE 
NETWORKS 

Our specific emphasis on the geometric limit of the 
macroscopic mechanical behavior of frictionless bead 
packings calls for an analysis of the geometry of sheared 
configurations, the first motivation of which is to explain 
the microscopic origins of macroscopic friction. Ulti- 
mately, a model should be sought which, unlike the anal- 
ogous one of Sec. IIIIE3[ would explicitly and quanti- 
tatively describe the mechanisms, involving instabilities 
and network rearrangements at the microscopic level, by 
which the material deforms and flows. Such goals will 
be only partly achieved here, since, leaving the detailed 
study of velocity correlations and strain mechanisms for 
future work, we focus on simple characterizations that 
are local in space and time. We also check here that the 
various microstructural variables studied, if measured in 
D-type simulations, approach their values observed in S- 
type ones in the limit of / ^ (at least in the large 
system limit). 

Packing geometry is classically described with a few 
state variables, among which the simplest ones are scalar: 
the volume fraction, the coordination number, as studied 
here in Section HV Al below. 

The much-studied distribution of contact force val- 
ues [4^ [47} is also determined in the present case 
(Section IIVB[) , and we check for effects of inertia and 
anisotropy. 

Under stress, or influenced by the history of their as- 
sembling process, the microstructure of grain packings 
develops anisotropic features, which are most often char- 
acterized with the fabric tensors, expressing statistics on 
orientations of normal directions at contacts, as stud- 
ied in Sec. IIV CI The critical state is microscopically 
characterized by stationary values of $, z, and fabric 
tensors, which are reached after a sufhciently large in- 
terval of monotonically growing strain in the quasistatic 
regime 142, ^,48,] • 



A. Coordination number 

The coordination number z strongly depends on / in 
steady state shear flow, and it is also affected by k. It 
decreases with increasing /, or with increasing k. As to 
the influence of C on z, it is notable for the largest / val- 
ues explored, but it decreases as the quasistatic limit is 
approached. Larger viscous damping coefficients increase 
the duration of contacts in shear flow, and thus produce 
slightly better coordinated networks on average. How- 
ever, the intensity of viscous forces becomes irrelevant in 
the quasistatic limit. According to our results, the de- 
pendence of z can safely be ignored for / < 10^"*. The / 
dependence of z is shown in Fig. 1121 

The coordination number of the equilibrated (S-type) 
anisotropic configurations is also very close to 6. This 
is a consequence of the isostaticity property of the 
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FIG. 12: (Color online) Coordination number 2; as a function 
of inertial number I, for k — ni (red square dots joined by 
dotted line, bottom data points) and k — K2 (blue crosses, 
top, dashed line). 



force-carrying structure (also called backbone [23|) of 
equilibrated sphere packings in the rigid limit - a re- 
markable property discussed in several recent publica- 
tions [m, [H, HI], which is specific to packings of rigid, 
frictionless and cohesionless spherical grains j4^, [s^ . 

Fig- HH shows that for quite low values of /, many con- 
tacts are lost (z is down to about 5 for / in the 10"'^ 
range) . The proportion po of rattlers increases with in- 
creasing /: Po is less than 1.5% in S simulations and for 
D simulations with / < 10~^ and k = ki, but is equal 
to 30% for / — 3.2 x 10^^ and k = ki. Our results are 
compatible with the theoretical value z = 6(1 — po) in 
the limit / ^ and k +00. Furthermore, it has been 
often observed that the grains only have a small number 
of contacts bearing large forces - this is the very reason 
why the "force chains" exist [U [5l|, [s^l • Consequently, 
as contacts carrying smaller forces are necessarily shorter 
lived, and tend to rarefy as / increases, the populations 
of grains with the largest local coordination are quickly 
depleted. 



B. Distribution of forces 

Fig. [13] is a plot of the probability distribution func- 
tion of the intergranular force normalized by the average 
force, for different values of /. The force distribution 
strongly depends on /: for / > 3.2 x 10^^ the probability 
distribution function p{f) (/ denoting the ratio of the 
normal force to the average value (-Fjv)) is monotonically 
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FIG. 13: (Color online) Probability distribution function of 
/ = N/{N) for I = 3.210"^ (red triangles), / = 3.210"^ 
(blue round dots) and / = 3.2 10"'' (black square dots). 



decreasing. For smaller values of /, p{f) has a maximum, 
around / = 0.5, and an approximately exponential de- 
cay for large values, as often observed in equilibrated 
granular packings [ll [H IH, [13, [HI, IH, H . The distri- 
butions obtained for the low values of / in D-type sim- 
ulations gradually approach the one obtained in S-type, 
equilibrated packings under maximum shear stress. The 
Kolmogorov-Smirnov test 55] can be used to detect the 
influence of parameters on the force distribution - the an- 
swer depending of course on the level of statistics of the 
available data. Based on 10 independent configurations 
of 4000 grains, it leads to the conclusions that no sig- 
nificant difference in force distribution could be detected 
between S-type results under maximum shear stress and 
D-type ones, and no influence of n either, provided the 
inertial parameter is small enough: / < 5 x 10~^, while 
some influence of C is only visible for /> 0.1. Our results 
are also compatible with a unique distribution, valid for 
maximum shear stress equilibrium configurations as well 
as for isotropic ones. 



C. Fabric 

Macroscopic friction is known to stem (at least par- 
tially) from the build-up of fabricjinisotropy in materials 
made of frictional beads or disks [4^[4^ . This connection 
is explored here with frictionless beads. 

Anisotropy of the tridimensional contact network can 
be characterized by the probability density function 
E{9, ip) of finding a contact with direction {6, ip). 6 is the 



colatitude angle and ip is the longitude angle of the spher- 
ical coordinates. E can be expanded in a series of spher- 
ical harmonics. The coefficients of the expansion are in 
one-to-one correspondence with the values of the fabric 
tensors, which are defined as the moments of the distribu- 
tion of normal unit vectors n on the unit sphere. Since a 
contact is left invariant by the parity symmetry f — r, 
E satisfies E{9, $) = E{'k - 0,ip + n). This means that 
the coefficients of odd order in the expansion in spherical 
harmonics are all equal to zero, and corresponds to the 
vanishing of all odd order fabric tensors. Coefficients can 
be computed from even order tensor products, viz. 



2k 



1 



2k 



cec 



(20) 



C denoting the set of Nc contacts, labeled with index 
c e C, where the normal unit vector is fic. 

Keeping only the lowest order of anisotropy, the expan- 
sion of E is restricted to the spherical harmonics of order 
two. Coefficients of the development are directly related 
to the value of the fabric tensor of order two, denoted by 



E{9,^) = 1/(471)+ Fi2d^yie,p) 
+ [Fii - F22)d,j.2_y2(e,(p) 

+ (F33- 1/3)42 (^^,^) + Fl3d.,(0,^) 

+ F23dyz{0, tf) + higher order terms (21) 

The constant 1 / (47r) corresponds to an isotropic distribu- 
tion and the next five terms of the development charac- 
terize the anisotropy of the material at the lowest order. 
Functions d are combinations of spherical harmonics of 
order two, with following expressions: 



dxy{0,p) 
dx2_y2{e,(p) 

dzAO,(p) 
dxz{0,(p) 
dvz{0,p) 



15 

8^ 
15 

16^ 
15 

16^ 

— sm o cos o cos ip 
47r 

15 • fl fl • 

— sm a cos a sm 09 

47r 



sin^ 9 sm(2ip) 
sin^ 9 cos{2ip) 
(3 cos^ 61- 1) 



(22) 
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Fabric tensor F is computed as a time average in the 
steady shear simulations. The numerical results show 
that F13 and F23 are always less than their respective 
statistical uncertainties, and can be considered as equal 
to zero, as requested by the symmetry in simple shear. 
We observe that F12 is always greater (by at least one 
order of magnitude) than the two other non negligible 
anisotropic coefficients, Fu — F22 and ^33 — 1/3. These 
two latter terms are below 2 x 10~^ for I < 10~^. Such 
low values are comparable with sample to sample fluctu- 
ations in equilibrated configurations. Thus, in the qua- 
sistatic limit, the anisotropy can be characterized by the 
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FIG. 14: (Color online) F12 as a function of inertial number 
I, with ( = 0.98 and N = 4000, for k = ki (red square 
dots connected by a dotted line) and for n — K2 (blue crosses 
connected by a dashed line). Both lines are power law fits of 
F12. 



sole F12 coefficient, the limit of which, as / 

-4 



0, 



evaluated as = —0.0165 ± 7. 10 ^ for k — ki and 



0.0156 ± 7. lO^'' for k = K2, with a fitting proce- 
dure. Like ^* and <&, F12 strongly varies with I (Fig.[T4|). 
and its dependence on / can be represented by a power 
law, with exponent ~ 0.36 for k = ki and ~ 0.37 for 

K = K2- 

F12 values measured in S-type equilibrated samples 
under maximum shear stress are influenced by system 
size N, very similarly to the static friction coefficient: 
larger values of |Fi2| are observed (typically ~ 0.02 for 
iV = 4000), but the excess over F12, the estimate from 
D-type simulations in the quasistatic limit, regresses as 
N increases, and the extrapolated macroscopic limit is 
compatible with the estimated values of F[2- This will 
be further examined in the more general context of the re- 
lationship between stress and anisotropics, for arbitrary 
stress tensors, in a forthcoming publication (56j . 

On changing ( from 0.98 down to 0.05, |Fi2| increases 
(correlatively with the decrease of z), by about 30% for 
/ ~ 10^^. This relative change is reduced to about 1% 
for / ~ 10^^ and the effect of C vanishes in the limit of 
/ ^ 0. 

Variations of F12 with parameters /, k and <^ are qual- 
itatively understood on noting that F12 is negatively 
correlated with the coordination number. If there are 
more contacting neighbors, on average, around a sphere, 
they are prevented by steric constraints from achieving 
highly anisotropic orientation distributions. This argu- 
ment, which with sirnple assumptions was made quanti- 
tative in 2D in Ref. [4^, thus explains that the increase 
of z observed as k is lowered tends to reduce |i^i2|, as 
observed on Fig. 1141 Similarly, the larger anisotropics 
observed away from the quasistatic limit are made pos- 



sible by the smaller number of contacts. The increase of 
IF12I with / is also due to the correlation of force intensi- 
ties with contact directions: on evaluating separately the 
fabric of the subnetworks corresponding to forces larger 
(or smaller) than the average contact force, one typically 
obtains, for / ^ 10~^, values of |_Fi2| twice as large (re- 
spectively: four times as small) as with the complete 
contact network. Contacts with small forces open if / is 
increased, and the remaining more strongly loaded ones 
are consequently more anisotropically oriented. 



V. DISCUSSION 

This work was devoted to the study of frictionless iden- 
tical spherical balls subjected to simple shear. The influ- 
ence of the three dimensionless quantities controlling the 
problem ~ inertial number /, stiffness number k and level 
of viscous damping - was carefully assessed and we ob- 
served that / has the most dramatic impact on the system 
behavior. Fluctuations of the measured quantities were 
shown to vanish for large systems. Consequently, the 
particular nature of the boundary conditions employed 
has no importance: for sufficiently large systems, fixed- 
volume simulations would lead to the very results we 
obtained with our stress driven numerical experiments. 
Particular attention was paid to the macroscopic geo- 
metric limit, that is the triple limit N +00, / — *■ 
and K +00. In this regime, the system behavior is 
governed by a succession of instabilities due to dynami- 
cal rearrangements of the contact network. A thorough 
investigation of such events remains an interesting, yet 
challenging, perspective. 

The existence of a nonzero macroscopic friction an- 
gle was evidenced by two different kinds of simulations - 
shear-rate controlled dynamic calculations (D-type simu- 
lations) and quasistatic stress-controlled calculations (S- 
type simulations). Whereas the dynamic friction angle 
is independent of the system size for N > 1000, the 
static friction angle (ps is very sensitive to the number 
of grains and is systematically greater than tpu for all 
studied sizes (TV < 8788) and (ps — Pd increases for 
decreasing N . This might be the reason why localiza- 
tion seems to occur more easily as the system size de- 
creases. In finite-size systems, the shear stress is a mul- 
tivaluated function of the strain rate in the quasistatic 
limit and the range of multivaluation increases with de- 
creasing N . Thus shear bands are more likely to appear 
in small systems [l^ 113, [El]. However, in the macro- 
scopic geometric limit, we found that both friction an- 
gles ips and (pD are equal within statistical uncertainties. 
In frictionless granular assemblies, all dissipation is due 
to viscous terms in contact forces, which therefore can 
be regarded 17[ as the physical origin of macroscopic 
friction. However, the value of the damping coefficient 
C, is irrelevant in the quasistatic limit since the amount 
of dissipated energy is geometrically determined. In the 
macroscopic geometric limit, we have seen that the shear 
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has no effect on the microscopic scalar quantities of the 
material (coordination number, distribution of forces), 
but it induces some structural anisotropy and a correla- 
tion between force intensities and contact orientations. 
We thus attribute the macroscopic friction angle to the 
shear-induced anisotropy of the material, as in the fric- 
tional case [s^. Ref. [5g| will show quantitatively that 
this is indeed the case. 

The result that ips — fD contrasts with observations 
on Lennard-Jones glasses at temperature T > [s^, HI] 
and on granular avalanches [60 . i61i] . Glass simulations 
show that the dynamic angle is less than the static one. 
This difference is linked to a stress overshoot visible on 
strain-stress curves. Similarly, in dense granular materi- 
als with friction, the shear stress goes through a maxi- 
mum before the steady state ("critical state") is reached, 
a feature which is absent in frictionless granular assem- 
blies (both states coincide in this case). Similar differ- 
ences {ips > fo) are reported for gran ular flows down 
inclined plane. Thus, in Refs. [goI. l6l|. 9stop{h) is less 
than 0start(^), where 9 is the inclination of the plane 
and h the thickness of the flowing layer in the station- 
ary state. The small thickness of the layer (typically less 
than ten grain diameters) and the intergranular friction 
are certainly responsible for this hysteresis. 

The stress-dilatancy interplay is a well known feature 
of granular materials. However, our simulations show 
that homogeneously sheared frictionless bead assemblies 
do not display any dilatancy in the macroscopic geomet- 
ric limit. In this limit, volume fraction $ remains equal to 
^'rcp during the whole time the material is sheared and 
the backbone stays isostatic in the rigid and quasistatic 
limits. This surprising lack of dilatancy can be intuitively 
understood in the light of the simple model presented in 
Sec. lIII iTSl We thus conclude that the steady state (crit- 
ical) volume fraction <l>c is equal to $rcp. 

The behavior of frictionless granular assemblies under 
arbitrary load directions will be the subject of a future 
work ^56] in order to gain a better knowledge of the yield 
surface and of the mechanical properties of such granu- 
lar systems under a small enough stress deviator (before 
failure). 

One motivation of the present work is the study of 
highly concentrated non-Brownian suspensions (Peclet 
number Pe = +oo), modeled as assemblies of ne arly 
touching grains bonded by a viscous lubricant [6^. [63l[6^ . 
Ideal lubrication effectively suppresses the tangent forces. 
Lubricated dynamics has already been employed as a 
means to obtain the force-carrying contact network of 
frictionless rigid particles, as the set of viscous bonds 
on which stresses concentrate [6^1 • Although crude, our 
current model should be able to reproduce the behav- 
ior of dense suspensions in the quasistatic limit. In this 
regime, the system evolves via a sequence of equilibrium 
states. At some point, the initial network is no longer 
able to sustain the imposed stress, it becomes unstable 
and a dynamic "crisis" occurs. Consequently, the evo- 
lution of the system is not quasistatic in the strictest 



sense (each point of the configuration space cannot be 
reached through a continuous series of equilibrium con- 
figurations). However, details of the dynamics are ex- 
pected to be irrelevant. Thus we expect that the same 
equilibrium states will be visited in the quasistatic limit 
by both frictionless granular systems and dense suspen- 
sions with frictional grains. According to the simple 
toy model of Sec. IIII E 3[ a dense suspension might be 
sketched by a slider moving on a bumpy surface in a 
media of viscosity rj. Close to the quasistatic limit, the 
most important parameter would be the dimensionless 
number rjj/ P. One may notice that it is very similar to 
the parameter /„ introduced by Cassar et al. that con- 
trols submarine avalanches in what they call the viscous 
regime j66| . Steady shear simulations evidenced that the 
material is still able to flow with a volume fraction ap- 
proximately equal to $* — <I>RCP — 0.64. This result is 
consistent with theoretical results pertaining to suspen- 
sions, where the volume fraction $* at which the viscos- 
ity of the suspension diverges is believed to tally with the 
random close packing volume fraction [6^ . However jit 
is not in agreement with the experiments exposed in [6|, 
where the value of <&* was found to be below 0.61. This 
discrepancy very likely originates in small scale features 
of the experimental system that are not accounted for a 
model of perfectly lubricated spherical beads. The be- 
havior of dense suspensions is known to be strongly im- 
pacted by short-range physics [68]. In the near future, we 
plan to study lubricated pastes with frictional contacts 
in the spirit of the simplified Stokesian dynamics scheme 
proposed by Ball and Melrose [13, [sl, [mI • 



APPENDIX A: CRYSTALLIZATION UNDER 
SHEAR 

Small samples, in both D and S-type simulations, tend 
to form strongly ordered structures under shear. This 
phenomenon, which do not occur for N > 1000, is briefly 
reported here. A more detailed study would be outside 
the scope of the present paper, and would require some 
investigation of the role of cell shape and boundary con- 
ditions, which is necessarily important in such small sys- 
tems. 

2 out of 3 S-type samples with N = 256 and stiffness 
level K2, and 2 out of 2 D-type samples with N — 500, 
/ — 3.2 X 10~* and k = ki present the following 
anomalies. First, solid fractions are considerably higher 
than 'I'rcp (and even more so considering the size ef- 
fect t21j, i23j on <i>) , with values approximatively equal to 
0.67 (see fourth line of Tab. IIIip . Apparent friction coef- 
ficients are also particularly large. A lower bound for the 
static macroscopic friction coefficient of S-type ordered 
samples is 0.4, whereas dynamic macroscopic friction co- 
efficient II* of D-type ordered samples for I = 3.2 x 10"'', 
K = Ki, and N — 500 may exceed by 20% the corre- 
sponding friction coefficient in bigger samples that do 
not experience any ordering. S-type samples also have 
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FIG. 15: Crystalline order induced by shear in one S sample 
with N ^256. 

very large coordination numbers, 8 < 2: < 9. This lat- 
ter characteristic is a clear indicator of partial crystalline 
order, as one necessarily has z < 6 in generically disor- 
dered situations. The denser crystal arrangements, face- 
centered cubic (fee) and hexagonal compact (hep) (and 
stacking variants thereof), have z = 12. For D simula- 
tions, anomalous values of $ and ai2 appear after strains 
of order 5. 



In order to detect crystalline order more quantita- 
tively, we use the standard order parameters Qe and 
Q4 employed in [13, [Hi, H [zll- Values of the pair 
((54, Qe) can be used to distinguish different local envi- 
ronments. In [2^, following [69J, the frequency of occur- 
rence of ranges of values (0.191 ± 0.05, 0.574 ± 0.05) and 
(0.097 ± 0.05, 0.485 ± 0.05), respectively corresponding 
to fcc-like and hcp-like configurations around one grain, 
were recorded. In the present case, most samples had 
very similar proportions of hcp-like and fcc-like local ar- 
rangements as in the RCP states studied in [2^: about 
12% of beads fall in the hep category, and fcc-like ones 
are virtually absent. The exceptions are the samples with 
anomalous, crystal-like properties, for which, while none 
of the beads has an fcc-like environment in that sense, 
the proportion of the hcp-like category raises to about 
60% in S samples and to 40% in D ones. 

A direct visualization. Fig. I15|, reveals strikingly or- 
dered configurations. A tentative conclusion to those 
preliminary observations is that the small samples tend 
to crystallize on somewhat shear-distorted hep lattices. 
One convenient characterization of order that is not sen- 
sitive to the distortion of crystalline patterns was sug- 
gested in [to], and used in '23|. With this method, more 
than 90% of the particles of the anomalous samples are 
declared to belong to crystalline regions. 
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